---
title: "Cloud seeding and rainfall"
author: "Christian Stratton and Meaghan Winder"
format:
  pdf:
    fontsize: 12pt
    fig-align: center
    mathspec: true
    geometry:
      - inner=1.5cm
      - outer=1.5cm
      - top=2.5cm
      - bottom=2.5cm
    header-includes: |
      \addtokomafont{disposition}{\rmfamily}
      \usepackage{fancyhdr, lastpage, framed, caption, xcolor, setspace, multirow, tabularray, framed, booktabs}
      \captionsetup[figure]{labelformat=empty}
      \pagestyle{fancyplain}
      \fancyhf{}
      \lhead{\fancyplain{}{STAT 0116: Introduction to Statistical Sciences}}
      \rhead{\fancyplain{}{Cloud seeding and rainfall}}
      \fancyfoot[RO, LE] {page \thepage\ of \pageref{LastPage}}
      \thispagestyle{plain}
      \usepackage{pdflscape}
      \fancypagestyle{mylandscape}{
        \fancyhf{} %Clears the header/footer
        \renewcommand{\headrulewidth}{0pt}% No header rule
        \renewcommand{\footrulewidth}{0pt}% No footer rule
      }
      \usepackage{stackengine,trimclip,scalerel}
      \savestack\eye{\rotatebox{90}{$^\circ\mkern-6mu\raisebox{1pt}{)}$}}
      \savestack\nose{\raisebox{3pt}{\scalebox{1}[-1]{\clipbox{0pt 1pt 0pt 0pt}{?}}}}
      \savestack\mouth{\rotatebox{90}{(}}
      \newcommand\Lenny{(\scalerel{\stackanchor[2pt]{\eye \nose \eye}{\mouth}}{)}}
editor: source
---

```{r setup, include = F}
# pull in packages and data
packs <- c('tidyverse')
lapply(packs, require, character.only = T)
options(show.signif.stars = F)
cloud <- Sleuth3::case0301
```

\doublespacing

## Introduction

Does massive injection of silver iodide into cumulus clouds lead to increased rainfall? Between 1968 and 1972, researchers in southern Florida conducted a randomized experiment in which clouds on 52 suitable days were randomly assigned to be seeded with silver iodide or left unseeded. Then, using radar, the researchers recorded the total volume of rainfall from each cloud. We seek to determine whether the average total volume of rainfall from the seeded clouds exceeds that from the unseeded clouds.

## Statistical procedures used

In this study, the response variable is the total volume of rainfall and the explanatory variable is whether the cloud was seeded; each of these variables were measured on the 52 clouds deemed suitable for cloud seeding, which represent the observational unit in this study. To address the research question, we conducted a one-sided two-sample t-test, which tests the following hypotheses:
$$
\begin{split}
H_0:& \mu_1 = \mu_2 \\
H_a:& \mu_1 > \mu_2
\end{split}
$$
where $\mu_1$ represents the true mean volume of rainfall in the seeded group and $\mu_2$ represents the true mean volume of rainfall in the unseeded group. This procedure relies on three assumptions: (1) independence of observations, (2) normality of errors, and (3) equality of variance between groups. In regards to independence, we do not see any egregious violations; since clouds were (presumably) only measured once, and observations were not collected sequentially over time, there is no reason to believe that rainfall across days would be related to one another. To examine the remaining assumptions, we plotted rainfall distributions for seeded and unseeded clouds (Figure 1), both on the raw scale and after a log transformation. The raw data show a strong right skew in both groups, with noticeably greater variability in the seeded group.

```{r eda, echo = F, fig.dim = c(6, 4), fig.align='center', fig.cap="Figure 1: Plot of rainfall volume by treatment on the log (bottom) and original (top) scales."}
# create dataframe for ggplot
rainfall <- c(cloud$Rainfall, log(cloud$Rainfall))
treatment <- rep(cloud$Treatment, 2)
resp_type <- factor(c(rep('original', nrow(cloud)), rep('logged', nrow(cloud))))

plot_df <- tibble(rainfall, treatment, resp_type)
plot_df$resp_type<- relevel(plot_df$resp_type, ref = "original")

# plot 
ggplot(plot_df, aes(x = treatment, y = rainfall)) +
  geom_violin() +
  geom_jitter(width = 0.1, alpha = 0.3, color = "lightblue") +
  facet_wrap(~ resp_type, scales = "free", nrow = 2) +
  coord_flip() +
  labs(title = 'Original and logged rainfall by seeding',
       x = 'Treatment',
       y = 'Observed rainfall (top) and Logged rainfall (bottom)') +
  theme_bw()
```


These features suggest violations of the normality and equal-variance assumptions. Additionally, the quantile-quantile plot of the raw data (see Appendix) deviates sharply from the hypothesized normal line, consistent with a violation of normality. Similarly, the larger spread in the seeded group, visible in both violin plots and the residuals vs. fitted plot (see Appendix), indicates heterogeneity of variance. Under both of these likely violations, our t-based inference is suspect at best.

To address this issue, we applied a natural log transformation to the rainfall data. This reduced skewness and produced violin plots with more symmetric shapes and comparable spreads. The QQ plot of the transformed data aligns well with the hypothesized normal line, and the residuals vs. fitted plot shows roughly constant variance (see Appendix). Thus, the transformed data are not inconsistent with the normality and equal-variance assumptions, and we proceeded with t-based procedures on the log-transformed values.

## Summary of Statistical Findings

```{r analysis, echo = F, results = 'hide'}
cloud$Treatment <- relevel(cloud$Treatment, ref = 'Unseeded') #relevel to match book
cloud.mod <- lm(log(Rainfall) ~ Treatment, cloud) #fit model
summary(cloud.mod) #summary
exp(1.1438) #estimate on original scale
exp(confint(cloud.mod)) #95% ci on original scale
```

We conducted a one-sided two sample t-test on 50 degrees of freedom, which yielded a test statistic of $t = 2.544$, with an associated p-value of $0.0141$. This provides strong evidence to suggest that the mean log-rainfall for the seeded clouds is greater than that of the unseeded clouds. On the original scale, it is estimated that the median total volume of rainfall for the seeded group is 3.14 times the median volume of rainfall for the unseeded group, with a 95\% confidence interval of 1.27 to 7.74.

## Scope of Inference

Since a random mechanism was used to determine whether clouds were seeded, we can conclude that the differences in rainfall that we observed were caused by the seeding. However, there is no reason to believe that the days on which clouds were seeded were randomly selected over the course of the four year experiment. Therefore, these results apply to only the 52 days deemed suitable to cloud seed, and only in the area in which the study was conducted.

## Acknowledgements and Author Contributions

Meaghan and Christian both wrote the report together. Meaghan made the exploratory graphics while Christian did the analysis with the help of ChatGPT 5.0, which was used to troubleshoot code to run the t-test.

\newpage

## Appendix

### Plots

```{r diag, echo = F, fig.dim = c(8,8)}
#create diagnostic plots for original and transformed data
par(mfrow = c(2,2))
plot(lm(Rainfall ~ Treatment, data = cloud), which = 1,
  main = 'Resid vs Fitted - Original Data')
plot(lm(Rainfall ~ Treatment, data = cloud), which = 2, 
  main = 'QQ Plot - Original Data')
plot(
  lm(log(Rainfall) ~ Treatment, data = cloud), which = 1, 
  main = 'Resid vs Fitted - Logged Data'
)
plot(
  lm(log(Rainfall) ~ Treatment, data = cloud), which = 2, 
  main = 'QQ Plot - Logged Data'
)
```

\singlespacing

### Code

```{r eda, eval = F, tidy = T}
```

```{r analysis, eval = F, tidy = T}
```

```{r diag, eval = F, tidy = T}
```